%% Calculate equilibrium


function [P_hat_n,P_hat_nj,P_hat_mnj,w_hat,X_hat_mnj,pi_mnjf1,pi_mnj1,...
    pi_c_mnj1,pi_c_nj1, pi_M1,pi_l1,pi_M_f1,pi_l_f1,PC_hat_n1,...
    pseudo_final_goods1,pseudo_intermediates1,Pi_hat_n]= ...
    EquilibriumCES_fun_approx_Profit(P_hat_mnj0,w_hat,pi_mnj0,pi_M0,pi_M_f0,...
    pi_l0,pi_l_f0,pi_mnjf0,pi_c_mnj0,X_hat_mnj0,pi_c_nj0,PC_n,zeta_mnj,X_mnj,Pi_hat_n0)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Produced for JdG,AL,IM by Christopher Evans at UPF  
%
% Program: Guess prices and wages, then iterate over prices, compare this 
% with the market clearing condition for the new guess of wages, then 
% iterate until convergence.
%
% Calls up the following functions:
%  1) PriceLoopCES_fun.m
% 	2) TradeShareCES_fun.m
% 	3) ExpenditureCES_fun.m
%	4) LabourMarketConditionCES_fun.m
% 
% Output: Converged expenditure hat, wage hat and prices 
% 
% NEW: For all code, we will change the subscripts to match those in
% master_notes3.pdf and the draft. I.e., we will change to {mn,ij} pairs,
% which represent country pair {mn} and sector pair {ij}. Global variables 
% changed and slight modification to file structure for calling up data.
% We also get rid of gamma_type variable. What previously _2.csv. files
% slightly as well as adjusting in our new directory.
%
% We also get rid of using the .(ISO{}) cell call up for the gammas and
% instead use a large matrix saved in Step2. This will help save on loops
% and speed up the proceedure.
%
% Last Updated: 5/02/2019 by IM & JDG
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global alphaT rhoT lambdaT etaT shockT rho sigma sigmaT lambda eta ...
    tol maxit D_hat_n Xi_hat_mnjf Xi_hat_mnj a_hat_f a_hat ...
    varphi sL_n sPi_n sD_n N J FI_J france oneminus_alpha_vector ...
    labour_share_PWT countrysecD firmsecD_sorted start start_sorted ...
    sum_by_country_dummy ISO alpha_WIOD_j alpha_WIOD_francej ...
    sector_algorithm 

w_max = 1; e = 1;

W = [w_hat'];
D = [];

fprintf('WAGE ITERATION=%d \n',e);

while (e <= maxit) && (w_max > tol)
   % Price Loop: Takes wage and price guess and loops until price
   % converges. Output is is the new P hat (change in prices) and wage hat
   % as well as expenditure hat (change in inputs/exports).
      
   [P_hat_mnj,b_hat,b_hat_f,P_M_hat,P_M_f_hat] = PriceLoopCES_fun_approx_Profit(P_hat_mnj0,w_hat,pi_mnj0,...
          pi_M0,pi_l0,pi_mnjf0,pi_M_f0,pi_l_f0);

   % Trade and labor Share: Creates the new trade (export and intermediate
   % input) and labor shares for the next period (t+1). Note that 
   % intermediate input and labor shares are not fed into P_hat_mnj 
   % function abovve, which is still a function of initial values only

  [pi_mnj1,pi_mnjf1,pi_c_mnj1,pi_M1,pi_l1,pi_M_f1,pi_l_f1] = ...
    TradeShareCES_fun_Profit(P_hat_mnj,b_hat,b_hat_f,P_M_f_hat,pi_M0,pi_l0,w_hat,pi_mnj0,...
    pi_M_f0,pi_l_f0,pi_mnjf0,pi_c_mnj0);

   % TESTING: reset firm-level pi's initial values (data) to see if we
   % converge
   
%    pi_l1=pi_l0;
%    pi_M1=pi_M0;
%    pi_l_f1=pi_l_f0;
%    pi_M_f1=pi_M_f0;
   
   % Expenditure: Creates a guess for the change in expenditure from t to 
   % t+1 and then loops over this until it converges. Follows the 
   % expenditure hat equation 23 from the master notes
   
   [P_hat_n,P_hat_nj,X_hat_mnj,pi_c_nj1,PC_hat_n1,...
       pseudo_intermediates1,pseudo_final_goods1,Pi_hat_n] = ...
      ExpenditureCES_fun_Profit(P_hat_mnj,pi_M1,pi_l1,w_hat,...
     pi_M_f1,pi_l_f1,pi_c_mnj0,PC_n,pi_c_nj0,...
   X_mnj,pi_c_mnj1,zeta_mnj,pi_mnjf1,pi_mnj1,X_hat_mnj0,Pi_hat_n0);
  
   % Labour Market Condition: Uses the labour market condition to update
   %                          the wage guess
      
   [w_hat_new] = LabourMarketConditionCES_fun_Profit(pi_l_f0, ...
            pi_l0,pi_l_f1,pi_l1,X_mnj,X_hat_mnj,...
            P_hat_n,pi_mnjf0,pi_mnjf1,pi_mnj0,pi_mnj1);  
 
        
   P_hat_mnj0 = P_hat_mnj;
   X_hat_mnj0 = X_hat_mnj;
   Pi_hat_n0 = Pi_hat_n;
   
   %w_max=sum(abs(w_hat_new-w_hat));
   
   w_diff=abs(w_hat_new-w_hat)./w_hat_new;
   w_max=norm(w_diff);

   
   if alphaT == 2
       
     % Convex iteration guess updater should be more stable and more likely
     % to converge as smaller jumps in the guess, but is slower to converge.
   
     w_hat=0.5*w_hat_new+0.5*w_hat;
     
   elseif alphaT == 4
       
     % Convex iteration guess updater should be more stable and more likely
     % to converge as smaller jumps in the guess, but is slower to converge.
   
    % w_hat=0.5*w_hat_new+0.5*w_hat;
      w_hat=w_hat_new; 
     % Updating the wage guess, this is the fastest but a convex iteration
     % updater may be more stable (smaller jumps of the guess)
   else
       
     w_hat=w_hat_new;
   
    
   end
   
%    % Update shares
%    pi_mnj0 = pi_mnj1;
%    pi_mnjf0 = pi_mnjf1;
%    pi_M0 = pi_M1;
%    pi_M_f0 = pi_M_f1;
%    pi_l0 = pi_l1;
%    pi_l_f0 = pi_l_f1;
%    pi_c_mnj0 = pi_c_mnj1;
    
   
%    % Small interation wage relative to previous step given a nu_W, where
%    % nu_W is the adjustment we might make
%    
%    nu_W = .5;
%    
%    w_diff = (w_hat_new-w_hat)./w_hat;
%    
%    w_hat_new = w_hat + nu_W*w_hat.*w_diff; % Iteration of new wage
%    
%    w_diff1 = (w_hat_new-w_hat)./w_hat;
%    
%    % Calculate distance as max between non-adjusted and adjusted new w_hat 
%    % and old w_hawt
% 
%    w_max = max(max(abs(w_diff)),max(abs(w_diff1)));
   
   % Set wage to new value to feed into next loop iteration
   
    W = [W w_hat_new']
    D = [D w_max]
   
%   w_hat = w_hat_new;
   
   e=e+1;

   fprintf('WAGE ITERATION=%d        Max wage distance=%d  \n',e-1,w_max);
end

end